transmission<-function(N=500,timesteps=501,mu=0.01,warmUp=300,top=NA,bias=0)
{
  agents <- 1:N
  traitCounter <- N + 1
  rawList <- vector("list",length= timesteps-warmUp)
  if (length(bias)==1){bias=rep(bias,timesteps)}
  for (t in 1:timesteps)
  { 	
    samplePool = agents
    if (length(unique(samplePool))>1)
    {
      sampleTraits = as.numeric(names(table(samplePool)))
      sampleProp = table(samplePool)/sum(table(samplePool))
      sampleTraitsProb = sampleProp^(1-bias[t]) / sum(sampleProp^(1-bias[t]))
      agents = sample(sampleTraits,size=N,replace=TRUE,prob=sampleTraitsProb)
    }
    index = which(runif(N)<=mu)
    if (length(index)>0)
    {
      newTraits<-traitCounter:c(traitCounter+length(index)-1)
      agents[index]=newTraits
      traitCounter=max(newTraits)+1
    }
    if (t>warmUp) {rawList[[t-warmUp]] = agents}	    
  }
  cases <- unique(unlist(rawList))
  rawMatrix <- t(sapply(rawList,instances,cases=cases))
  obs.div <- 1 - apply(rawMatrix,1,function(x){sum(c(x/sum(x))^2)})
  if (is.na(top)) {top = min(apply(rawMatrix,1,function(x){sum(x>0)}))} 
  zMatrix=turnover(rawMatrix,top=top)
  z = apply(zMatrix,2,mean)
  z.frame=data.frame(y=1:top,obs.z=z)
  aN = 1.38 * (mu^0.55) * N^0.13 
  z.frame$exp.z.Bentley=z.frame$y*sqrt(mu)
  z.frame$exp.z.EvansGiometto= aN*(z.frame$y)^0.86	/ 2
  z.frame$zfitN= z.frame$obs.z/aN
  modN=lm(log(zfitN+0.000001)~log(y),data=z.frame)
  b=as.numeric(coefficients(modN)[2])
  
  return(list(x = b, obs.div = obs.div))
  
  #remove variables and do garbage collection
  rm(list = c("modN", "z.frame", "aN", "z", "zMatrix", "rawMatrix", "cases", "t", "rawList", "warmUp", "agents", "traitCounter", "b", "obs.div",
              "newTraits", "index", "N", "mu", "sampleTraitsProb", "sampleTraits", "sampleProp", "samplePool", "timesteps", "bias"))
}

instances <- function(x,cases)
{
  x <- c(x,cases)
  return(table(x) - 1)
}

turnover <- function(mat,top)
{
  z<-matrix(NA,nrow=c(nrow(mat)-1),ncol=top)
  for (i in 1:c(nrow(mat)-1))
  {
    t1 = names(sort(mat[i,],decreasing=TRUE))
    t2 = names(sort(mat[i+1,],decreasing=TRUE))
    z[i,]=sapply(1:top,function(x,t1,t2){return(x-length(intersect(t1[1:x],t2[1:x])))},t1=t1,t2=t2)
  }
  return(z)
}

reScale = function(x,lo,hi)
{
  return(((hi-lo)/(max(x)-min(x)))*(x-max(x))+hi)
}